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Practical flux prescriptions for gamma-ray burst 
afterglows, from early to late times 
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ABSTRACT 

We present analytic flux prescriptions for broadband spectra of self-absorbed and 
optically thin synchrotron radiation from gamma-ray burst afterglows, based on one- 
dimensional relativistic hydrodynamic simulations. By treating the evolution of critical 
spectrum parameters as a power-law break between the ultrarclativistic and non- 
relativistic asymptotic solutions, we generalize the prescriptions to any observer time. 
Our aim is to provide a set of formulas that constitutes a useful tool for accurate fitting 
of model-parameters to observational data, regardless of the dynamical phase of the 
outflow. The applicability range is not confined to gamma-ray burst afterglows, but 
includes all spherical outflows (also jets before the jet-break) that produce synchrotron 
radiation as they adiabatically decelerate in a cold, power-law medium. We test the 
accuracy of the prescriptions and show that numerical evidence suggests that typical 
relative errors in the derivation of physical quantities are about 10 per cent. A software 
implementation of the presented flux prescriptions combined with a fitting code is 
freely available on request and on-line f. Together they can be used in order to directly 
fit model parameters to data. 

Key words: hydrodynamics radiation mechanisms: non-thermal radiative transfer 
shock waves gamma-ray burst: general. 



1 INTRODUCTION 

Gamma-ray bursts (GRBs) are believed to be produced by 
powerful relativistic outflows resulting from the catastrophic 
death of massive stars (Woosley 1993), or the merger of 
two compact objects (Eichler et al. 1989). The burst itself 
(prompt emission) likely arises from internal shocks occur- 
ring due to the variability of the central engine (Rees & 
Meszaros 1994; Sari & Piran 1997), while the afterglow emis- 
sion comes from the interaction of the same outflow with the 
medium surrounding the burster (Rees & Meszaros 1992; 
Paczyhski & Rhoads 1993). Although the dominant radia- 
tion process behind the prompt emission is not yet clear, 
it is well established that the afterglow radiation is domi- 
nated by synchrotron emission from shock-accelerated elec- 
trons (Meszaros & Rees 1993; van Paradijs et al. 2000). 

The prompt emission is typically very brief and concen- 
trated at high energies. On the other hand, afterglows are 
often visible over many more orders of magnitude both in 
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time- and frequency-space (see Meszaros 2006 for an exten- 
sive review of GRB research). Thus, studying the afterglow 
radiation allows us to put a multitude of constrains both 
on the microphysics (e.g. the fraction of internal energy go- 
ing to the magnetic fields and the power-law accelerated 
electrons) governing the shocked plasma (Spitkovsky 2008; 
Sironi & Spitkovsky 2009), as well as on the basic physical 
parameters describing the phenomenon macroscopically, like 
blast-wave energy, density and structure of the surrounding 
medium. 

It is these macroscopic parameters that determine the 
dynamical evolution of the outflow. However, a full analytic 
description of the dynamics is only possible when the spa- 
tial component of the four-velocity of the outflow /?7 is either 
much greater (Blandford & McKee 1976) or much smaller 
(Sedov 1959) than 1. Therefore, relativistic hydrodynamic 
(RHD) simulations (Kobayashi et al. 1999; Meliani et al. 
2007; Zhang & MacFadyen 2009; De Colle et al. 2012) are 
the most accurate means of studying the intermediate dy- 
namical regime linking the ultrarelativistic and Newtonian 
solutions (see however Huang et al. 1999). Van Eerten et al. 
(2010a) have numerically studied the lightcurves of outflows 
advancing through all three dynamical regimes and have 
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shown that the transition is slow, i.e. deviations from the 
expected relativistic behaviour appear well before the New- 
tonian asymptotes are reached, mainly due to the changing 
adiabatic index of the shocked gas. For typical burst param- 
eters (isotropic blast-wave energy Ei BO — 10 52 ergs, ambi- 
ent medium number density no = 1 cm~ 3 ) the Sedov- Taylor 
scalings set in at a few thousand days, observer time, imply- 
ing that an appreciable portion of the afterglow (typically 
around hundreds of days) emanates from outflows with dy- 
namics that cannot be described analytically by either of 
the two asymptotic solutions. 

Soon after the discovery of the first afterglows (Costa 
et al. 1997; Groot et al. 1997) efforts were made to calcu- 
late broadband synchrotron spectra and light curves as a 
function of burst parameters (Wijers et al. 1997; Sari et al. 
1998; Panaitescu & Kumar 2000). The common way to do 
this is by tying the dynamical evolution of the blast-wave 
in regimes where this is feasible to radiation models that, 
according to the jump conditions at the shock front, cal- 
culate the resulting spectra. Despite the success of early 
efforts in capturing general features of the observed spec- 
tra, the progressive refinement of the used models has led 
to very different estimates of the physical parameters of in- 
dividual bursts. For example Wijers & Galama (1999) and 
Granot & Sari (2002) have both fitted GRB 970508 and 
their derived values differ up to 3 orders of magnitude. Fur- 
thermore, the applicability of most of these models is re- 
stricted to a particular dynamical phase and only recently 
have there been a few attempts at addressing the entire evo- 
lution of spectra and light curves through the performance 
of simulations (Zhang & MacFadyen 2009; van Eerten et al. 
2012; Wygoda et al. 2011; De Colle et al. 2012). Even so, 
these models do not contain a treatment of self-absorption 
(apart from van Eerten et al. 2012), necessary to model low- 
frequency observations with e.g. the Expanded Very Large 
Array EVLA (Perley et al. 2011), the Low- Frequency Ar- 
ray LOFAR (Morganti et al. 2011) and the upcoming Karoo 
Array Telescope MeerKAT (Booth et al. 2009) and Square 
Kilometre Array SKA (Carilli & Rawlings 2004), and do 
not provide flux-prescriptions. Van Eerten et al. (2012) do 
provide a broadband fit code, but it requires the use of a 
parallel computer network. 

The purpose of this work is to provide accurate analytic 
flux-prescriptions, based on one-dimensional RHD simula- 
tions, that are applicable to both the ultrarelativistic and 
Newtonian phase but also, and perhaps more importantly, 
to observer times when the outflow is transitioning from the 
former to the latter. Apart from the typical, initially ultra- 
relativistic outflows of GRBs, the formulas we present are 
applicable to Newtonian as well as relativistic (Soderberg 
et al. 2010) outflows from supernova explosions in the adi- 
abatic phase (Chevalier 1977, 1982; Draine & McKee 1993) 
and mildly relativistic outflows originating from binary neu- 
tron star (NS) mergers, expected to produce detectable elec- 
tromagnetic (EM) counterparts to gravitational wave detec- 
tions (Nakar & Piran 2011; Metzger & Berger 2012). They 
can also be applied to relativistic outflows resulting from 
the tidal disruption of stars by a super-massive black hole 
(Bloom et al. 2011; Metzger et al. 2012), under the limiting 
assumption of quasi-spherical outflow. The presented model 
naturally accounts for the exact shape of the synchrotron 
spectrum (including self-absorption, but ignoring cooling) 



and the structure of the blast-wave. Furthermore, it can be 
applied to a range of power-law density structures of the cir- 
cumburst medium, a possibility previously studied by van 
Eerten & Wijers (2009) and De Colle et al. (2012). This al- 
lows for modelling of more complex environments, expected 
on a theoretical basis (Ramirez-Ruiz et al. 2005) and de- 
duced observationally (Curran et al. 2009). With such a tool 
a light curve can be fitted without the need of costly simu- 
lations and the restrictions of models specialising in specific 
dynamical phases, or preset structures of the circumburst 
medium. In order to obtain the flux-prescriptions we com- 
bine three elements: (1) analytic formulas for flux-scalings 
during the Blandford-McKee and the Sedov- Taylor phases, 

(2) one-dimensional, hydrodynamic simulations, using the 
adaptive mesh refinement code AMRVAC (Meliani et al. 2007; 
Keppens et al. 2012), that span the whole range of the dy- 
namics (from ultrarelativistic to Newtonian velocities) and 

(3) a radiative-transfer code that uses simulation snapshots 
and a parametrisation of the microphysics to calculate in- 
stantaneous spectra. 

This paper is organised as follows: in section 2 we briefly 
describe the setup of the performed simulations and the sub- 
sequent calculations of spectra and light curves. In section 
3 we present formulas that describe the flux as a function 
of physical parameters in both the relativistic and the New- 
tonian phase of the outflow. That includes specifying the 
flux at any given power-law segment, as well as a descrip- 
tion of the sharpness of the spectral breaks that occur at 
critical frequencies. We then proceed in section 4 to connect 
the two dynamical regimes (relativistic and Newtonian) by 
treating the transition from the former to the latter as a 
prolonged temporal break the characteristics of which can 
be linked to the physical parameters of the burst and its 
environment. In section 5 we describe how one can make 
use of the flux-prescriptions to obtain spectra at any given 
time. We also show comparisons between spectra based on 
simulations and spectra constructed using the provided pre- 
scriptions. Finally, we present an application of this model to 
mildly relativistic outflows from binary neutron star mergers 
in order to assess the recent predictions of Nakar & Piran 
(2011) concerning the detectability of the produced radio 
signals. In section 6 we discuss our results and the implica- 
tions of this work for GRB afterglow models. 



2 NUMERICAL TREATMENT 
2.1 Simulations 

We have made use of the AMRVAC adaptive- mesh-refinement 
numerical code to run a series of simulations for different val- 
ues of physical parameters. These simulations span a wide 
range of the four- velocity at the shock front, from ultrarel- 
ativistic values (~ 70) down to ~ 0.05. 

In total 7 different simulation runs were used to arrive 
at the presented prescriptions. They can be characterized 
by the blast-wave energy E52 (in units of 10 52 erg), start- 
ing Lorentz factor of the shock Tin, maximum radius of the 
simulation-box i? max , slope of the power-law density distri- 
bution of the surrounding medium k and value of the number 
density at 10 17 cm, no. In Table 1 we present the values of 
these parameters for each run. 
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Table 1. Parameters of simulations used to derive the 
flux prescriptions. Left column enumerates the per- 
formed simulations. 



Sim 


E52 


k 


n 


r in 


-Rmax (cm) 


1 


1.0 


0.0 


1.0 


60 


3 • 10 19 


2 


0.04 


0.0 


4.0 


60 


10 19 


3 


1.0 


0.5 


1.0 


60 


4 • 10 19 


4 


0.4 


0.75 


0.5 


28 


5 • 10 19 


5 


1.0 


1.0 


1.0 


60 


4 • 10 19 


6 


1.0 


2.0 


1.0 


70 


10 19 


7 


0.01 


2.0 


1.0 


10 


7- 10 19 



2.1.1 Resolution 

The ID simulation box typically extends from 10 16 cm up 
to a few times 10 cm, although both limits were modified 
accordingly for physical models with different external den- 
sity profiles and blast- wave energies. Utilizing the adaptive- 
mesh-refinement approach in AMRVAC we have used a max- 
imum of 20 refinement levels which set the effective reso- 
lution of the grid. For 120 cells at the lowest refinement 
level, this amounts to a resolution of ~ 4.77 ■ 10 11 cm per 
cell. For comparison, the corresponding width of the ini- 
tial BM shell for Simulation 1 is (van Eerten et al. 2010b) 

■Rin/(6r2 n ) =5 . 1Q 12 

cm, with 7?i n denoting the radius of the 

shock. 

We have checked for convergence against runs of differ- 
ent refinement levels (see also van Eerten et al. 2011, De 
Colle et al. 2012) and have also checked our results against 
theoretically predicted values in the early part of the out- 
flow when the resolution demands are the highest. They have 
been found in good agreement. 



2.2 Radiative-transfer code 

The snapshots generated by the simulation runs were post- 
processed using a radiative-transfer code (van Eerten & 
Wijers 2009; van Eerten et al. 2010a). During the post- 
processing an array of beams is created and propagated at 
the speed of light through the three-dimensional generalisa- 
tion of the ID snapshots towards the observer. The elements 
of this array take the value of the specific intensity I v . At 
each step the solution to the equation of radiative transfer 
is applied to each beam and the value for the intensity is 
updated through the equation 

I v =I e- T " +S„(l-e- T »), (2) 

where Io is the value of the intensity at the previous step, t„ 
is the optical depth and S v the source function. We note that 
the optical depth of a single simulation-cell can be larger 
than unity. 

The array of intensities is constructed so that the po- 
sitions of its elements lie on the surface from which light- 
signals arrive at the observer simultaneously. The density 
of the beams is determined through an adaptive-mesh ap- 
proach ensuring sufficient resolution. Once all beams have 
crossed the entire blast- wave integration of the intensity over 
the surface yields the flux 

A 

where d is the luminosity distance to the observer, (1 + z) 
the cosmological correction and A the surface defined by the 
beams. 

In the case of on-axis jets and spherical outflows, like 
the ones considered in this study, eq. (3) can be reduced to 
a ID integral due to the axisymmetry of the source (Granot 
et al. 1999). Our approach makes use of that symmetry and 
a ID integral is solved numerically to calculate the observed 
flux. 



2.1.2 Equation of state 

In all the simulations we have used a 'realistic' equation 
of state (EOS) with an effective adiabatic index (Meliani 
et al. 2004) that lies between the ultrarelativistic and non- 
relativistic limits (4/3 and 5/3, respectively): 

rad.cff = § - 1 (1 - . (i) 

In the above equation p'c 2 is the comoving rest-mass en- 
ergy density and is weighed against the total energy density 
(including rest-mass) of the gas u' . 

This Synge-type (Synge 1957) EOS has also been used 
in van Eerten et al. (2010a), where the effects have been 
analysed and comparisons to constant r at j, ff have been 
made. In short, its effect on the observed flux is that of 
a very gradual transition from values close to (but not at) 
those corresponding to an ultrarelativistic EOS to values 
approaching those of a non-relativistic EOS. 



3 FLUX-PRESCRIPTIONS IN THE 

ASYMPTOTIC DYNAMICAL REGIMES 

In this section we demonstrate how to combine blast-wave 
dynamics in each of the self-similar regimes with synchrotron 
radiation theory to arive at scalings that describe the ob- 
served flux as a function of frequency and time. 



3.1 Shock dynamics 

We first outline the dynamics of the outflow. As mentioned 
in section 1, one can obtain power-law scalings for the four- 
velocity as well as the mass and energy densities right behind 
the shock front as a function of blast- wave radius and time 
in the two extreme regimes of dynamical behaviour of the 
afterglow. 
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3.1.1 Ultrarelativistic phase 

In the ultrarelativistic (also known as Blandford-McKec, 
hereafter BM) phase these scalings take the form (Bland- 
ford & McKee 1976): 



72 CC T sh oc t~ 



ci oc F sh m, 



n 2 oc T sh ni . 



(4) 
(5) 
(6) 



where 72 is the Lorentz factor of the shocked plasma mea- 
sured in the lab frame (which coincides with the frame 
of the surrounding medium), describes its internal en- 
ergy density and ri2 its number density. Here and through- 
out this paper the quantities e and n will be measured in 
the comoving frame of the fluid they are describing. F sh 
is the Lorentz factor of the propagating shock wave, while 
ni(r) — no (r/ro) is the density of the unshocked medium 
surrounding the burster as a function of radius. In all the 
results presented the characteristic distance ro is put at 
10 17 cm. The time H' appearing in eq. (4) is the lab-frame 
time and is to be distinguished from the observed arrival 
time of light signals which is affected by light travel-time 
effects. While in the BM phase we can assume r ~ ct for 
the radius of the shock. 



3.1.2 Newtonian phase 

At the phase where the outflow has become Newtonian (also 
known as Sedov- Taylor phase, hereafter ST) a similar ap- 
proach can be taken to describe its kinetic and thermody- 
namic evolution. Dimensional analysis implies 



r(t) oc 



Et 2 
m(r) 



1/5 



(7) 



for the scaling of the radius of the shock as a function of 
time, blast-wave energy and structure of the surrounding 
medium. This leads to 



/3(t)<xt 



e2 oc 1 5 ~ 



712 OC Til, 



(8) 
(9) 
(10) 



dr 



where /3 = — — is the bulk velocity of the shock, in units of 
cat 

c. 



3.2 Optically thin and self-absorbed synchrotron 
radiation 

The next step is to combine these scalings with formulas that 
calculate optically thin as well as self-absorbed synchrotron 
radiation. We assume that electrons are the primary radiat- 
ing particles and express their post-shock energy distribu- 
tion as a power-law N(E) oc E~ p . The lower limit of the 



distribution E m = 7, 
chrotron frequency v' m — 



c 2 corresponds to a comoving syn- 



"7n 



Qc B' 



since, where Q e and 



4-7T rn c c 

ra e are the electron charge and mass, respectively, B' is the 
comoving value of the magnetic field and a is the pitch angle 
between magnetic field and velocity of the electron. In the 
case of optically thin radiation the flux at a given frequency 
v' , in the comoving frame will be 



F' v , oc (p-l)NB'Q(y m ), 



(11) 



where N is the total number of power-law accelerated elec- 
v 1 

trons, y m — — and 



v' 



i-p r 

= X 2 / 
Jo 



y 2 V(y)dy. 



(12) 



The function Q contains all the spectral information. The 
function V appearing in eq. (12) is the synchrotron function 
F(x) (Rybicki & Lightman 1986) integrated over all pitch 
angles, for an isotropic pitch-angle distribution. 

In the case of optically thick synchrotron radiation the 
comoving flux is given by 



F' 



(13) 



where j' v and a' v are the comoving emissivity and absorption 
coefficient, respectively, and r 2 is a measure of the radiating 
surface. The expressions for j' v and a' v have the form 



3u oc (p - l)£n 2 B' Q(y m ), 



(14) 



e? 1 B' v'- 2 Q{p+l,y m ), (15) 



, P-1 2 p + 2 f2 2 - 
p-2 

where ^ is the fractional number of electrons accelerated 
to a power-law distribution, e c is the fraction of internal 
energy carried by the accelerated electrons and Q(p + 1, 2/ m ) 
is evaluated using eq. (12) by replacing p with (p + 1). The 
effect of absorption is the introduction of another critical 
frequency in the spectrum, v & . The ordering of f m and v & 
determines the shape of the spectrum. In Fig. 1 and 2 the 
two different spectra are represented schematically in order 
to illustrate the break frequencies as well as the slopes of 
the power laws that they connect. 

To arrive at the relations above we have demanded that 
the distribution of the electrons obeys 



f 

J E, 



N(E) dE = £n 2 



and 



/ N(E)EdE = e e e 2 , 

J 



(16) 



(17) 



where we have used the implicit condition that p > 2. 

In this study we ignore the effect of electron-cooling 
on the spectra since the fast timescales associated with it 
translate to distances much shorter than the typical size of 
a simulation-cell. This work focuses on observer times when 
the influence of cooling on the observed spectra is negligible. 
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Figure 1. Spectrum 1. Normalised form of the spectrum when 



surrounding medium, respectively. The former quantity (p) 
determines the electron distribution, everywhere behind the 
shock front, for a given set of thermodynamic parameters, 
while the latter (k) affects the structure of the decelerating 
blast-wave. Two standard values for k are often assumed 
in the literature, namely (constant density medium) and 
2 (constant stellar wind environment). However, fits to k 
(Yost et al. 2003; Curran et al. 2009) often indicate differ- 
ent conditions, motivating us to use it as a free parameter. 
Values for the factors of the calibrating polynomial are then 
derived by demanding that they satisfy the system of equa- 
tions resulting from runs of models with different p and k. 
The range of values we have explored are [2.1, 3] for p and 
[0, 2] for k. Therefore this is also the range under which the 
presented prescriptions are applicable. 

Having put all of the ingredients together, the equation 
describing the flux at any given power-law segment of the 
synchrotron spectrum, either in the BM or the ST phase, 
has the general form 



10° 




v a 


„ 10 2 
-o 

CD 
CO 




r ^\ 


(normal 
o 


v m/ 


v (i-p>v 


^ 10* 


r 




10~ 8 








10 6 


10 8 10 10 10 12 10 14 
v(Hz) 


Figure 2 


Spectrum 2. 


Normalised form of the spectrum wl 



F v = Q.O. h(p) £«« e *« 4* ng» EM tls (l+*)"<5?> 



< fa- 



3.3 General form of fiux-scalings 

Equations (4)-(10) allow us to calculate the conditions right 
behind the shock front as a function of time. From these 
equations we can compute instantaneous spectra by utilizing 
standard formulas for synchrotron radiation. However, eq. 
(4)-(10) do not specify the structure of the shocked plasma 
well behind the shock. Such a specification would have al- 
lowed us to convolve the different parts of the outflow that 
contribute to the observed radiation at any given observer 
time. Such an approach has been taken, for example, by 
Granot & Sari (2002). 

Our approach is based on the fact that in the self-similar 
regimes the scaling-behaviour of the emitted radiation can 
be calculated by considering a homogeneous slab that obeys 
the scalings of the shocked fluid right behind the shock. How- 
ever, in order to correctly calibrate the scalings (i.e. provide 
the correct flux-levels) one has to capture the shock struc- 
ture behind the front and the most reliable way to do this 
is by simulations. 

The calibration is done by introducing a polynomial in 
terms of p and k - the spectral index of the power-law accel- 
erated electrons and the index describing the structure of the 



(18) 



where 



logC^ i = go + g p p + g pp p 2 + g k k + g kk k 2 (19) 

and h(p) is a function of p, different for each power-law seg- 
ment. It takes the following values: 



Hp) = ^ 



(p - 2) 3p + 2 



(p-l)(p + 2) 3p-l' 

1 G(p) 
(p + 2) G(p+1)' 



p — 1 I p — 2 



where 



3p — 1 \p — 1 
(P - 1) G( P ) 



G(p) = 



-2/3 



p-2 



p-i 



F v oc v 
F v oc 

F v oc 
F v oc 



(20) 



rtf + !)r(! + i§)rt! 



r(| + f) (p + i) 



(21) 



G(p) as well as other factors appearing in eq. (20) originate 
from the limiting behaviour of Q(x). For details see van 
Eerten & Wijers (2009). 

Equation (18) shows all the possible physical depen- 
dencies of the flux that this model is taking into account. 
We have introduced the fraction of internal energy carried 
by the magnetic field 6b- This quantity, along with e e , £ 
and p constitute a group describing the microphysics of the 
shocked electrons and enter via synchrotron theory. There- 
fore, the exponents cjj, q e , qa and the prefactor h(p) remain 
the same regardless of the dynamics of the outflow. On the 
other hand there are two quantities describing the burster 
and its environment: no and the blast- wave energy E52 (mea- 
sured in units of 10 52 erg), two quantities describing the fre- 
quency (i^obs) and the time (t bs) of the observation and two 
more describing the cosmological distances usually associ- 
ated with GRBs: the redshift 2 and the luminosity distance 
d,28 (in units of 10 28 cm). 

We note that the inclusion of £ in our description of the 
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microphysics has the implication that one cannot uniquely 
determine the values of all model parameters at once. This is 
a consequence of the degeneracy of the used model which for 
a set of primed parameters E' 52 = E 52 /f, n' =no/f,e' e = 
/e e , £b = / £ b, £' = /£, produces the same spectrum as 
the set of unprimed parameters. This degeneracy was first 
pointed out by Eichler & Waxman (2005) and can also be 
seen in eq. (24), (25) and Tables 9-12 presented in Section 
4. As a result, a value for one of the parameters must be 
assumed during fitting in order to determine the others. 

All the g-exponents appearing in eq. (18) are deter- 
mined analytically. They are in general unique for a par- 
ticular power-law segment in a given dynamical phase of 
the outflow. This also holds for all the g-factors appearing 
in eq. (19). Their values, however, are determined by match- 
ing them to numerical results (for a variety of (p, k) values) 
and solving the resulting systems of equations. They are in 
fact the calibration of the flux-scalings. 

3.4 Flux-scalings 

3.4-1 Flux-scalings during the BM phase 

A very similar approach has been taken by van Eerten & 
Wijers (2009). These authors have explored optically thin 
synchrotron radiation from relativistic outflows, taking into 
account all the possible spectra that result from either fast 
or slow cooling (Sari et al. 1998). Here we expand on that 
by including self-absorption. Table 2 contains the values of 
the g-exponents (analytically derived dependencies), while 
Table 3 contains the values of the (/-factors (numerically 
determined calibration). 

The values of the (/-exponents are in agreement with 
the formulas presented in Granot & Sari (2002), apart from 
the fact that we have chosen to include an extra parameter 
£. The normalisation of the flux-scalings results in slightly 
lower fluxes (of order 30%) compared to Granot & Sari 
(2002) a difference which can be attributed to the varying 
adiabatic index of the simulations (van Eerten et al. 2010a). 
Below u m we have ignored stimulated emission associated 
with a population inversion of the electron distribution at 
the low-energy limit, something which Granot &: Sari (2002) 
have included in their model. That accounts for a factor of 
approximately [3(p + 2)/4] difference in flux between those 
predictions and the present ones. 

3.4-2 Flux-scalings during the ST phase 

For the Newtonian phase of the outflow, we repeat the same 
procedure as in the BM phase. The analytically derived de- 
pendencies are presented in Table 4 while the factors of the 
calibrating polynomials are presented in Table 5. We note 
that the values of q v and q t are in agreement with those 
presented in Frail et al. (2000), apart from their equation 
(A18) where the scaling for ^ obs <S v m is in error. This error 
also appears in van Eerten et al. (2010a). 

3.5 The sharpness of spectral breaks 

In practice the spectral breaks are not infinitely sharp as 
shown in Fig. 1 and 2 but show a gradual transition from 
one power-law index to another. To complete our description 



Table 2. (/-exponents in the BM phase. Analytically derived 
g-exponents for self-absorbed and optically thin synchrotron 
radiation in the BM phase. The quantities on the left corre- 
spond to the different physical parameters on which the flux 
depends (see eq. 18). Each column describes their values for 
a given power-law segment. 
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Table 3. (/-factors in the BM phase. Numerically deter- 
mined (/-factors for self-absorbed and optically thin syn- 
chrotron radiation in the BM phase. The quantities on the 
left correspond to different factors of the normalizing poly- 
nomial (see eq. 19). Each column describes their values for 
a given power-law segment. 
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Table 4. (/-exponents in the ST phase. Analytically derived q- 
cxponents for self-absorbed and optically thin synchrotron radi- 
ation in the ST phase. The quantities on the left correspond to 
the different physical parameters on which the flux depends (see 
eq. 18). Each column describes their values for a given power-law 
segment. 
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Table 5. g- factors in the ST phase. Numerically deter- 
mined g-factors for self-absorbed and optically thin syn- 
chrotron radiation in the ST phase. The quantities on the 
left correspond to different factors of the normalizing poly- 
nomial (see eq. 19). Each column describes their values for 
a given power-law segment. 
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Table 6. s-factors in the BM phase. Numerically deter- 
mined s-factors (see eq. 23) for all possible breaks in the 
BM phase. Each column describes a specific break, with 
the two associated spectral indices denoted on top. 
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of instantaneous spectra, we need to provide a formula for 
the sharpness of spectral breaks. 

An approach commonly used (Granot & Sari 2002; van 
Eerten & Wijers 2009) is to describe the flux close to a break 
by the following equation (Beuermann et al. 1999) 



10- 



Simulation-based spectrum 
Fitting function 

s=3.0 




v(Hz) 

Figure 3. The effect of sharpness on the flux close to a spectral 
break. This fragment of the simulation-based spectrum focuses 
on the flux around u m . The best fit is shown along with two more 
curves that have the same parameters but different sharpness. 
Simulation-based spectrum has the following model parameters: 
E 52 = 1, n = 1, p = 2.5, k = 0, £ = lO" 2 , e c = 10" 1 , e B = 
10~ 2 , d 2s = 1, z = 0.56, t obs = 100 days. 



Table 7. s-factors in the ST phase. Numerically deter- 
mined s-factors (see eq. 23) for all possible breaks in the 
ST phase. Each column describes a specific break, with the 
two associated spectral indices denoted on top. 
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F v (i> oba ) = A 



-1/s 



(22) 



where (vq,A) are the coordinates of the meeting point of 
the two power laws associated with the break, a\ and a 2 are 
the asymptotic power-law indices before and after the break, 
respectively, and s is the so called 'sharpness parameter'. We 
have performed x 2 - mm hnization fitting in logarithmic space 
to obtain values of s for specific runs and used those to arrive 
at a description of the sharpness in terms of a polynomial 
of p and k. This polynomial has the general form 



s = so + s p p + Sk k + Skk k 



(23) 



Its factors have been determined by solving the system of 
equations resulting from the application of eq. (22) to models 
with different p and k parameters. In Tables 6 and 7 we 
present the values of so,s p ,Sk and Skk in the BM and the 
ST phase, respectively. 

In Fig. 3 a best fit to the shape of the spectrum around 
v m is shown, for one of the run models. We also plot the flux 
for two different values of s to illustrate the notable effect it 
can have on flux levels. 



4 THE TRANSRELATIVISTIC REGIME 

The results of the previous section constitute a full descrip- 
tion of the possible synchrotron spectra (ignoring cooling) 
during the ultrarelativistic and non-relativistic dynamical 
phases of the afterglow evolution. However, we have not 
yet addressed a large portion of the afterglow's overall be- 
haviour, namely the transrelativistic regime. During this 
stage the dynamics deviates considerably from the BM so- 
lution, without having settled yet into the ST solution. As 
mentioned in section 1, this phase of the afterglow typically 
spans a few orders of magnitude in observer time, while there 
is no full description of its dynamics, even in the simple, 
spherical case. 

An approach we have investigated and found useful is 
that of treating the transrelativistic phase as a 'break' dur- 
ing which the temporal evolution of the spectrum's critical 
parameters displays a smooth transition from the relativis- 
tic power-law behaviour to the non-relativistic one. These 
parameters could be, for example, the values of the flux at 
every possible power-law segment. However, the same level 
of accuracy can be achieved by using the positions of the 
critical frequencies and the flux at one of them, instead. 
Based on values for these parameters one can construct the 
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Table 8. g-factors for -FJn_BM and 
-fm-ST- Numerically determined g- 
factors for Fm (measured at i/ m i), 
both in the BM and ST phase. 
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spectrum (Sari et al. 1998; Wijers & Galama 1999) because 
the slopes of the power-law segments are known for a given 
ordering of v & and f m . 



Table 9. / n and g-exponcnts for critical frequencies in the BM 
phase, while u a < v m (spectrum 1). The g-exponents carry the 
analytically derived dependencies, while the / n -factors carry the 
flux-calibrating C* po i and h(p) (see eq. 18, 19 and 20). 
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4.1 Peak flux 

A convenient frequency to measure the flux is at v m of spec- 
trum 1. This is because we can assume that the bulk of 
the electrons radiate most of their power at that frequency. 
Using eq. (11) in combination with the scalings for the dy- 
namics in each of the two extreme phases of the outflow, we 
find for the flux at v m 



1 4 8-3fc -k s _ k 

F m . BM = Cpoi^I £™ t 2 a it k) {l+z)*P=>* d^, 

(24) 



in the BM phase and 



i?n-ST = Cp i £ e£ n 



7 8-3fc 3-2fc 

2(5 -fc) ,-,2(5-fc) j. 5-fc 



*«*T a+z)—d^, (25) 



in the ST phase. The (/-factors of C po i for both dynamical 
phases are presented in Table 8. 



Table 10. / n and g-exponents for critical frequencies in the 
BM phase, while f m < u a (spectrum 2). 
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4.2 Critical frequencies 

The behaviour of the critical frequencies can easily be de- 
duced (in either the BM or ST phase) by equating the flux 
formulas on both sides of a spectral break. For each of Vm 
and f a there will be two such expressions corresponding to 
the two possible spectra for the two different orderings of 
the frequencies. In general, the value of a critical frequency 
will be given by the following formula 



fn 



ft 



E 



52 obs obs 



(26) 



The numerical factors / n result from equating the fluxes of 
the power-laws at each side of the spectral break. Tables 9 
and 10 summarize the formulas for the critical frequencies 
in the BM phase for the two different possible spectra. For 
the ST phase we repeat the same procedure and summarize 
our results in Tables 11 and 12. 

For clearer presentation we have labelled as v^i and z/ m i 
the critical frequencies v & and i/ m , respectively, when v a < 
v m (i.e. when spectrum 1 applies), while they are labelled as 
^ a 2 and v m 2 in the opposite case when spectrum 2 applies. 



Table 11. / n and g-exponents for critical frequencies in the ST 
phase, while < v m (spectrum 1). The expressions contained 
in f n need to be evaluated using the formulas applicable to the 
ST regime. 
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Table 12. / n and g-exponents for critical frequencies in the ST 
phase, while i/ m < v a (spectrum 2). 
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4.3 Evolution of critical parameters 

A practical way of describing the temporal evolution of the 
parameters needed to construct a spectrum at any point is 
that of a smoothly broken power law. We can use eq. (22), 
this time characterizing a temporal break in the following 
manner 



$(lobs) 



^obs 

lo~ 



^obs 



-1/st 



(27) 



where (to, A) is the meeting point of the asymptotes and $ 
is the value of any of the critical parameters. In this version 
of eq. (22) ai and ai are the BM and ST slopes, respectively. 
We can rewrite the above equation in the following way 



.) = (*i5+<"£? 



-1/s 



(28) 



where both $bm and $st have to be evaluated at t ohs . 

In the previous sections we have established not only 
the scalings of the critical parameters we wish to follow, 
but their actual values as a function of observer time in 
both extreme dynamical regimes. This allows us to insert 
them directly into eq. (28), where the only unknown left is 
the sharpness s t . As in the case of spectral breaks we have 
performed ^-minimization fitting in logarithmic space and 
have arrived at a description of 's t ' in terms of a polynomial 
of the following form 



St = so + s-p p + Spp p 2 + Sk k + Skk k 2 



(29) 



Results for the values of these s t -factors are presented in 
Table 13. 

An example fit of F m is shown in Fig. 4. The behaviour 
of F m displays clear deviations from the BM scalings already 
before 100 days observer time, for Ei BO = 10 52 erg, no = 1 
and k = 0. It settles to values sufficiently close (within 10%) 
to the ST solution at around 5000 days. The duration of 
the transrelativistic regime is represented in s t , for a given 
set of physical parameters. Table 13 demonstrates that the 
sharpness of every parameter is generally unique. Based on 
that we conclude that duration and features of the tran- 



Simulation-based values of F m 
Fitting function 
BM-slope 
ST-slope 




10"' 



10' 



10' 
Us ( da y s ) 
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10' 



Figure 4. A broken power-law fit to the evolution of F m . Plotted 
are the BM and ST asymptotes. Simulation-based data have the 
following model parameters: £52 = 1, no = 1, p = 2.5, k = 0, £ = 
1(T 2 , e c = lO -1 , e B = 1CT 2 , (2 2 8 = 1, 2 = 0.56. 



Table 13. s t -factors for the evolution of critical parameters. 
Numerically determined st-factors (see eq. 29) describing 
the evolution of critical parameters (break-frequencies and 
maximum flux) from the BM to the ST phase. Each column 
describes a specific parameter, the name of which is denoted 
on top. 
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"The s-factors for 1/^2 have not been determined by 
solving the system of equations resulting from measur- 
ing its value in different models but comprise a rather 
heuristic approach that minimizes the deviations from 
the numerically determined values. 



srelativistic phase will be manifested differently across the 
spectrum. 

4.4 Evolution of the sharpness of spectral breaks 

Our findings so far enable us to determine the values of the 
critical frequencies and F m at any given time for a burst of 
given physical properties. This allows for an accurate cal- 
culation of the flux at any given power-law segment of the 
spectrum. What is left to specify is the flux close to a spec- 
tral break for a general t ohB . To achieve that we need to 
provide a quantitative description for the evolution of the 
sharpness of spectral breaks from the BM to the ST phase. 

As it turns out the sharpness of every spectral break 
follows a 'characteristic path' as it evolves from the relativis- 
tic values to the Newtonian ones. This path is qualitatively 
independent of the physical properties of the burst and is 
unique for every spectral break. In Fig. 5 we present these 
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log/ 



NR 



Figure 5. Evolution of sharpness for all possible spectral breaks. 
Before t; and after if, s resumes its standard BM and ST values. 
The sharpness of the break around u m i shows the most complex 
pattern, by declining initially during the transrclativistic regime 
and then rising again to meet its ST asymptote. The value of s m i 
at tNR is ~ 0.18 smaller than the BM value. 



paths for all possible breaks. From the data gathered we 
have identified three timescales that are represented in this 
figure: U, 4nr and if. In fact, we can simplify things further 
by setting t\ — tNR/ 100 and tf — 10 tNR which is generally 
valid up to a few percent, regardless of the physical param- 
eters of a burst. 

Determination of £nr carries (as in the case of fluxes 
and critical frequencies) an analytic and a numerical com- 
ponent. The analytic part is motivated by considerations of 
the dynamics and is similar to other estimates of an observer 
time marking the transition to the Newtonian phase (Livio 
& Waxman 2000; Piran 2004). Specifically it is identified as 
the observer time at which the shock Lorentz factor drops 
to the value of 2, following the BM solution. 



£nr 



fc-4 
l3-k 



where Anr 



17 



8 7rm,r 



(30) 



Including the numerical calibration the expression for 
tNR takes the form 



t NR = 10 1 



fc-4 
}3=fc" 

r~ fc^ NR 



3-fc 



(1 + z) days. 



(31) 



Like all other timescales describing a transition from 
BM to ST, t N R scales as (-E 5 2/n ) 1/(3 ~ fc) (1 + z), as we ex- 
pect from dimensional analysis (van Eerten & MacFadyen 
2012). The same holds for to appearing in eq. (27) for all 
critical spectrum parameters. However, the actual value of 
to for every parameter is influenced by the flux-calibrating 
polynomials Cp \ (eq. 19) and is therefore in general unique. 

The equations, tables and plots of sections 3 and 4 carry 
all the information necessary to construct a spectrum at any 
given time, based on given values for the relevant physical 
parameters that we have discussed in this work. In the fol- 
lowing section we demonstrate how these results can be used 



to calculate the observed flux at any given frequency and 
time. 



5 USING THE PRESCRIPTIONS 

In this section we focus on the practical side of this work 
which is to construct spectra at any given observer time 
based on values for the physical parameters of a burst. These 
parameters we repeat here for clarity: £52 (isotropic blast- 



wave energy in units of 10 erg 



no 



(number density at 



10 cm), p (index of the electron power-law distribution), 
k (index of the density distribution of the matter surround- 
ing the burster), £ (fraction of accelerated electrons), e c and 
£b (fractions of internal energy assigned to the relativistic 
electrons and magnetic field, respectively) , v h s and t b s (fre- 
quency and time of observation) , d,2s and z (luminosity dis- 
tance in units of 10 28 cm and redshift). 

The task of constructing a spectrum out of the pre- 
sented formulas can be divided in four parts. The first is to 
obtain the values of the critical frequencies and the max- 
imum flux at a given observer time. The next step is to 
determine the shape of the spectrum and its general char- 
acteristics, i.e. values of the flux away from the breaks, for 
each of the two possible spectra. The third step is to assign 
the appropriate sharpness parameters to all spectral breaks. 
The fourth and final step is to use eq. (32) and (37) in order 
to calculate the observed flux at a given frequency. 

Each of the two spectra is described by a single equa- 
tion, at a given observer time. This equation should es- 
sentially represent a mathematical formulation of a double- 
broken power law. One of the ways to achieve that (Granot 
& Sari 2002) is to use a heuristic formula that combines eq. 
(22) with a factor that assigns a second break at a different 
frequency. Then the whole spectrum can be described by the 
following expression 



F„{v ohs ) = A 



+ 



fobs 
V0 



-1/s 



1 + 



h(a 2 -03) 



-l/h 



(32) 



The first two terms on the right-hand side of the above equa- 
tion describe the first break as usual, while the third term 
describes the second break. We have introduced v\ the fre- 
quency of the second break, 03 the slope of the third power- 
law segment and h the sharpness of the second break. 

Equation (32) is exact only when the two breaks of the 
spectrum are sufficiently away from each other so that the 
power law connecting them is apparent, even for a small 
range of frequencies. When this is not the case, this equation 
provides an approximation to the real spectrum (Granot & 
Sari 2002). 

In Fig. 6 we present a flowchart of the basic steps to- 
wards creating a spectrum of the emitted radiation at any 
given observer time. The details of each step can be found 
in the following subsections of the text. 
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Calculate F m , f m , i/ a (Section 5.1) 



Work out the critical parameters for 
each type of spectrum (Section 5.2) 



Calculate the appropriate sharpness 
for each spectral break (Section 5.3) 



Use eq. (32) and (37) to cal- 
culate the flux (Section 5.2) 



Figure 6. Flowchart showing the basic steps of the presented 
method. References are given to specific sections of the paper 
where the steps are described in more detail. 



5.1 Values of F m , v m and i/ a 

In order to decide which of the two possible synchrotron 
spectra is valid (see Fig. 1 and 2) one needs to calculate 
all i; a i, "a2, "mi and v m 2- While all these frequencies can 
in principle be calculated independently, we opt for a dif- 
ferent method. Requiring that the flux at both ends of the 
spectrum (in the power- laws v 1 and u 1 ^ 1 ^ 2 ) be the same 
regardless of the type of spectrum, one can show that only 
three of the four frequencies are independent. The relation 
between them is 



"a2 



10 3p-l 1 

,.3(4+p) 3(4+p) 4+p 
-al "ml "m2 ■ 



(33) 



We have expressed in terms of the others because 
its transrelativistic profile is the one deviating more from 
the broken power-law approach. This way inconsistencies 
that may arise during a spectral transition, due to the over- 
specification of the evolution of the spectrum, are avoided 
and the flux in the leftmost and rightmost power laws is 
independent of the ordering of the critical frequencies. 

Equation (28) should be used to calculate i/ a i, "mi 
and Vm2 for a given set of physical parameters and for the 
same observer time, before applying eq. (33) to obtain i/ a 2- 
$BM(tobs) and $sT(£ bs) assume the corresponding forms of 
the critical frequencies, as those are expressed in eq. (26) 
and Tables 9, 10, 11 and 12. The value of s t is given by eq. 
(29) and the relevant entries of Table 13. 

Equation (28) should also be used to calculate F m at the 
same observer time as the critical frequencies. The asymp- 
totic expressions are presented in eq. (24) and (25) and are 
to be evaluated using eq. (19) and Table 8. Having found 
the values of all 5 critical parameters at the same observer 
time, we can now start constructing the spectrum. 



5.2 Shape and flux-normalisation of the spectrum 

In the case of spectrum 1 the characteristic synchrotron fre- 
quency f m i lies on the optically thin part of the spectrum. 
Consequently, its value is affected by radiation from the 
whole blast-wave. In the case of spectrum 2 self-absorption 
allows only for the front to contribute to the flux close to 
f m 2 • In practice the values of v m i and v>m2 are always close 
to each other (within the same order of magnitude). There- 
fore, we conclude that the location of i/ m on the spectrum is 
mostly determined by the conditions at the front and is only 
slightly affected by what the optical depth of the blast-wave 
is at that frequency. 

On the other hand, values of ^ a i and ^ a 2 can differ sub- 
stantially with respect to each other. However, in most cases 
they will both be either smaller or bigger than z/ m i and 
" m 2- We will be referring to these cases as definite order- 
ing, whereas all other cases will be referred to as indefinite 
ordering. When z/ a i , ^ a 2 < "mi , " m 2 the spectrum will have 
the form of Fig. 1 (spectrum 1), while if z/ m i, " m 2 < "ai, "a2 
that of Fig. 2 (spectrum 2). In the case of indefinite order- 
ing, the actual positions of the two critical frequencies on 
the spectrum are very close to each other signaling a spec- 
tral transition, typically from spectrum 1 to spectrum 2. In 
terms of observer time, the time-span of this transition is 
relatively small. Let t tr be the observer time when indefinite 
ordering sets in. The duration of this transtion (time-span 
of indefinite ordering) is typically a fraction of t t r- During 
that time the choice of critical frequencies affects the flux 
across the spectrum by factors of order unity. 

In practice, it is preferable to always take the values 
suggested by both spectra into account, through a consis- 
tent weighing method. This way glitches that may appear in 
the produced lightcurves when switching from one spectrum 
to another are avoided. Instead, the lightcurves' behaviour 
smoothly progresses from the early-time spectrum 1 con- 
figuration to the late-time spectrum 2. The weight of each 
spectrum is represented by a power-law dependence in time 
that contains a characteristic timescale tm P related to the 
observer time at which the spectrum is transitioning from 
spectrum 1 to spectrum 2. This characteristic timescale can 
be estimated numerically by solving eq. (27) for both frri 
(the time at which z/ m i = "ai) and £t2 (the time at which 
" m 2 = "a2), and defining 



teip = /flip • max(tTi, te)- 



(34) 



We have found that the value of /flip that results in smaller 
deviations across the parameter space of [p, k] is 1.6. 
The weights of spectrum 1 and 2 can be written as 



Wi = 



W 2 



(^obs/iflip) 1 



(iobs/iflip) 1 + (tobsAflip) 
(tobs/tflip) 



(iobs/iflip) 1 + (tobs/iflip) 

The flux at a given frequency will then be 



(35) 
(36) 



log F = Wi • log Fi + W 2 ■ log F 2 , (37) 
where Fi and F 2 are the fluxes calculated at that frequency 
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through spectrum 1 and 2 (see subsections 5.2.1 and 5.2.2), 
respectively. 

5.2.1 Spectrum 1: u a < v m 

This is the asymptotic case where both z/ a i and f a 2 are 
smaller than z/ m i and ^ m 2. Consequently the positions of 
the critical frequencies are given by v & = i/ a i and v m — v m \. 
The parameters of eq. (32) get the following values: 

(i) Vq = Vax 

(ii) fl = fml 

(iii) A — F m (Vza) 1/3 

(iv) ai = 2 

(v) a 2 = 1/3 

(vi) o s = (l-p)/2 

5. ,2. 2 Spectrum 2: v m < u a 

For spectrum 2 in the asymptotic limit both i/ a i and ^ a 2 are 
bigger than v m \ and ^ m 2. Thus, the positions of the critical 
frequencies are given by i/ m = u m 2 and ^ a = i/ a2 . However, in 
spectrum 2 the flux at v m is not F m ; that would be the case 
if it were not for absorption. We can use this fact to first 
calculate the flux at v & . Although in the present spectrum- 
configuration the actual position of u m is given by i/ m 2, it 
is I/mi that we should use for obtaining the flux at v & . The 
variables become: 

(i) Vq = U m2 

(ii) V! = i/ a2 

/ \ (l-p)/2 / x 2.5 

(-) ^ = te) te) 

(iv) ai = 2 

(v) a2 = 2.5 

(vi) a 3 = (l-p)/2 

5.3 Sharpness parameters of spectral breaks 

The only two parameters left to specify in eq. (32) are s 
and h, the sharpness of the first and the second break in 
the spectrum, respectively. In order to assign the proper 
sharpness to each break we first have to compare tabs to ^nr 
(see eq. 31). There are four distinct cases, as illustrated in 
Fig. 5: 

5.3. 1 t Q }) S <c ti 

In this case all sharpness parameters attain their BM values 
as these are given in Table 6. 

5.3.2 t obs >t f 

In this case all sharpness parameters attain their ST values 
as these are given in Table 7. 

5.3.3 ti < t a l,s < tNR 

In this case breaks v m 2 and ^ a2 retain their BM sharpness. 
The other two exhibit some evolution towards the corre- 
sponding ST values. Namely, the sharpness around f m i will 
be 



Smi = 0.09 logf ) + Si , (38) 

V *obs / 

where Si is the sharpness of the particular break in the BM 
regime. 

The sharpness around v & \ will be 

Sai = - — jj — -log! — 1+S;, (39) 

where s\ and Sf are the sharpness parameters at the BM and 
ST phases, respectively. 

5.3-4 tNR < t t s < tf 

In this final case all breaks exhibit a sharpness evolving to- 
wards its ST value. For z^ m i the sharpness will be given by 

Ami = (Sf - Si + 0.18) l°g(^) + s u (40) 

while for u & \ the value of the sharpness is still given by eq. 
(39). 

The sharpness around ^ a 2 will be given by 

«b2 = (*-*) log +Si, (41) 

while for f m 2 we find a similar result 

s m2 = (s f - Si )log(|^+ S i. (42) 
5.4 Examples of results 

We have described a practical implementation of our re- 
sults to construct spectra at any given time, based on val- 
ues for the physical quantities characterising the burst. We 
now show comparisons between simulation-generated spec- 
tra and spectra that have been constructed using the pro- 
vided flux prescriptions. 

In Fig. 7 a comparison between a simulation-based spec- 
trum and an analytic one is shown. In all power-law seg- 
ments and the linking breaks, the flux-prediction is never 
more than 10% off compared to the simulation-based data. 
We can translate these deviations into relative errors for the 
values of the physical parameters. This we do by adjust- 
ing their values so that those deviations vanish in particular 
regimes of the spectrum. For the blast-wave energy (-B52) 
the error ranges from 3% up to 15% depending on which 
power-law segment (or spectral break) one uses for the com- 
parison. In the case of no the maximum error is 15%. For p 
the difference is of order 0.05, while for k it is of order 0.08. 
For e c and eb the error is of the order 10% while for £ it 
reaches up to 30%. 

In Fig. 8 we present another comparison between a 
simulation-based spectrum and a constructed one. This one 
was chosen for exhibiting one of the largest deviations we 
have encountered. While the self-absorbed part of the spec- 
trum is matched well by the constructed spectrum, flux in 
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Figure 7. A typically good match between an analytically con- 
structed spectrum and one based on a simulation. Both are taken 
at 100 days, v & lies at ~ 10 7 Hz and u m at ~ 10 12 Hz. Model 
parameters for both spectra are: £52 = 1, no = 1, p = 2.3, k = 
0, £ = 10- 2 , e c = 10-\ e B = 10~ 2 , d 2& = 1, z = 0.56. 



5.5 Application to mildly relativistic outflows 

Recently Nakar & Piran (2011) have discussed the radio- 
signal following the ejection of spherical, Newtonian or 
mildly relativistic outflows expected from binary neutron 
star mergers. They estimate that due to the low initial 
Lorentz factors of these outflows, their deceleration (and 
entry to the ST phase) will be manifested at idec ~ 60 days 
observer time, for £52 = 0.01, no = 1, k = 0, $ ~ 1 (initial 
velocity). This is also the time at which optically thin emis- 
sion at u bs = 5 GHz will peak in the range 0.01 — 0.1 mjy, 
for a distance of the source in the range 1 — 3 Gpc. Here, we 
test these estimates using the prescriptions presented in this 
paper. 

The model we have developed in this study is based on 
(and therefore, applicable to) outflows that are initially ul- 
trarelativistic. Thus, it is not obvious that it can be used 
to model non-relativistic outflows. Order of magnitude cal- 
culations in the lab frame can illustrate the limitations. A 
relativistic outflow of coasting Lorentz factor Ti will slow 
down after sweeping mass Fj times smaller than the mass of 
the ejecta (Rees & Meszaros 1992). This will happen at a 
time 



Constructed spectrum 
Simulation-based spectrum 




Figure 8. An example of a constructed spectrum that shows 
relatively large deviations from the numerical result in the opti- 
cally thin part of the spectrum. Both spectra are taken at 500 
days. Um lies at ~ 1 Hz and v& at ~ 10 7 Hz. Model parameters 



are: £52 = 1, no 
10~ 2 , d 2S = 1, 2 = 



0.56. 



2.5, k = 0.5, ? = 1, e e 



10- 



the v 1 p,/2 segment differs by ~ 25%. The corresponding er- 
rors in the derivation of values for physical quantities are the 
following: for E 52 up to 16%, for n 10% - 90% (flux in the 
gment depends very weekly on no for these model 
parameters), for £ up to 35%, for e e 15%, for eb 30%, while 
for p and k we find differences up to 0.1 and 0.2, respectively. 

We stress that these deviations are not with respect 
to a best-fit value but are indicative of how much every 
parameter should be tweaked to match fluxes in individual 
power-law segments of the spectrum. More often than not 
such a tweak would actually produce a rather bad fit overall. 
Thus the deviations we have listed may be viewed as an 
upper limit to what a broadband fit would produce. 



\4-7rpiC 5 



(43) 



From that point onwards the outflow will decelerate accord- 
ing to the BM solution (r oc t~ 3 ' 2 ) becoming Newtonian 
(r ~ 1) at 



£n = 



/ 3E 



\4npic 5 
The corresponding radius is 

rN = £n c. 



1/3 



(44) 



(45) 



Equations (44) and (45) effectively mark the onset of the ST 
phase. 

In the case of sub- and mildly relativistic outflows the 
deceleration time (also marking the transition to the ST 
phase) occurs when the swept mass is comparable to the 
rest mass of the ejecta 



^dec 



3E 



2irpic 5 

At tdec the shock is at a radius 



1/3 



Pi 



■5/3 



3i tdec C. 



(46) 



(47) 



From eq. (44)- (47) it is clear that as /3i — > 1 the onset of 
the ST phase for Newtonian outflows approaches that of the 
relativistic analog with the same energy. This implies that 
fast (v ~ c) outflows (regardless of the Lorentz factor) have 
no memory of their history from idec onwards. Therefore we 
can apply the ST scalings of the flux-presciptions to a mildly 
relativistic outflow (as the one considered by Nakar & Piran 
2011) at observer times t bs ^ tdec- In the sub-relativistic 
case (Pi <§C 1) eq. (46) and (47) imply that the outflow will 
decelerate later and at a greater radius. Nevertheless, the 
ST scalings of the flux-prescriptions apply at observer times 
£obs 3> tdec, i.e. sufficiently later than the deceleration time. 

For the application to mildly relativistic outflows we 
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Figure 9. Synchrotron spectra of a mildly relativistic outflow at 
distances of 1 and 3 Gpc, taken at 60 days observer time. The 
grey shaded area represents the 5 GHz band of the EVLA with a 
bandwidth of 3.5 GHz and a 4cr detection threshold of ~ 10 fijy 
for lhr integration (Perley et al. 2011). The value of p for both 
spectra is 2.5. The other physical parameters have the following 
values: £52 = 0.01, no = 1, k = 0, £ = 1, e e = ee = 0.1. 

have set £52 = 0.01, no = 1, k = for the macroscopic 
parameters of the blast-wave and its environment and £ = 
1, £e = £b = 0.1, for the microphysics, while v h s = 5 GHz 
and t bs = 60 days. The electron spectral index p is varied 
within the range 2.1 — 3.0, while for the distance we have 
taken the two extreme values of 1 and 3 Gpc. 

In accordance with Nakar & Piran (2011) we find that 
i^obs is in the optically thin part of the spectrum for all cases. 
In this regime the flux increases monotonically for an in- 
creasing p and for p = 3.0 it is about four times higher than 
the p = 2.1 case. For d 28 = 0.31 (~ 1 Gpc) the flux at 5 GHz 
lies in the range 0.015 — 0.06 mjy, depending on the value of 
p. At d'28 = 0.93 (~ 3 Gpc) we find the flux to be always be- 
low 0.01 mjy, albeit marginally for relatively high values of 
p. This is illustrated in Fig. 9 where two spectra are shown 
corresponding to distances of 1 and 3 Gpc. They are both 
taken at t obs = 60 days, for a characteristic value of p = 2.5. 

Binary neutron star mergers are believed to be the 
progenitors of short GRBs (see Nakar 2007 and references 
therein) and are observed in a variety of environments, like 
elliptical, spiral and irregular galaxies (Berger 2009). A con- 
siderable fraction of them, however, appear to be host-less, 
occuring in the intergalactic medium and thus surrounded 
by a much more tenuous gas than that commonly found 
inside galaxies (Berger 2010). We have repeated the calcu- 
lation of the spectrum from a spherical outflow resulting 
from a NS-NS merger for a surrounding medium of density 
no = 10 -3 cm~ 3 , where the lower density results in a later 
onset of the ST phase at ~ 400 days. Keeping all other pa- 
rameters constant we find that the radio signal at 5 GHz will 
be detectable by the EVLA up to a distance of ~ 100 Mpc. 

The implication of these results is that moderately en- 
ergetic outflows (E — 10 50 erg) expected to accompany NS- 
NS mergers (Rezzolla et al. 2011) can produce synchrotron 
radiation detectable by the EVLA from distances up to 
~ 1 Gpc, larger than the detection horizon of the upcom- 
ing versions of gravitational-wave detectors (Nakar & Piran 
2011). This assumes that the density of the matter surround- 



ing the merger is of the order 1 cm -3 . These radio signals will 
peak at timescales of the order of a few months, if the cor- 
responding outflows have initial velocities close to the speed 
of light. In the case of more tenuous circumburst media the 
ST timescale grows and the detection horizon of the EM sig- 
nal drops accordingly. The presented flux-prescriptions are 
applicable throughout the ST phase of these outflows. 



6 DISCUSSION 

We present analytic flux-prescriptions, for broadband syn- 
chrotron spectra originating from GRB outflows, suitable 
for fast and detailed modelling of the afterglow phase. 
They are applicable throughout the evolution of observed 
afterglows, during which external shocks are the domi- 
nant source of particle acceleration, and account for the 
exact shape of the synchrotron spectrum, including self- 
absorption, but ignoring cooling. These prescriptions are 
based on high-resolution, one-dimensional, hydrodynamic 
simulations performed using the adaptive mesh refinement 
code AMRVAC. To obtain spectra we have employed a ra- 
diation code that solves the equation of radiative transfer 
through the evolving blast-wave as this is determined by 
the simulations. The presented formulas carry two compo- 
nents. The first is derived analytically and expresses the 
dependence of the flux on relevant physical parameters 
(E52, n , p, k, £, e e , tB,^obs, t obs , d 2 8, z), while the second 
component reflects the calibration that the results of simu- 
lations have introduced to the flux-levels. 

For each asymptotic dynamical regime (BM and ST) we 
provide prescriptions for the flux at every power- law segment 
but also at frequencies close to spectral breaks. These are 
modelled as smoothly broken power-laws with the sharpness 
of the break given in terms of the structure of the surround- 
ing medium (k) and the electron distribution (p). During 
the transrelativistic regime we find that the values of critical 
frequencies (i/ m , z/ a ) and peak flux (-F m ) of the synchrotron 
spectrum show a gradual transition from the asymptotic 
power-law behaviour in the BM phase to the correspond- 
ing one in the ST phase. This fact has allowed us to model 
their temporal profiles as smoothly broken power-laws. For 
every parameter (F m , v m and i/ a ) we provide formulas de- 
scribing the sharpness of these breaks in terms of p and k. 
In order to model the evolution of a spectral break's sharp- 
ness, we have recognized the unique pattern that each break 
exhibits. We have introduced £nr whose derivation is based 
on considerations of the outflow dynamics. The result is a 
set of analytic expressions that extend the applicability of 
the flux-prescriptions to any given observer time. 

An element of this study worth emphasising is the in- 
clusion of k (representing the structure of the circumburst 
medium) as a fitting parameter. This is motivated by the 
fact that environments of stars with variable mass-loss rates 
(such as massive stars, prime candidates for long GRB pro- 
genitors) can have structures more complex than the usually 
assumed k = or 2 (Ramirez-Ruiz et al. 2005) and fits us- 
ing k as a free parameter do not exclude such a possibility 
(Yost et al. 2003; Curran et al. 2009). As can be seen in 
Tables 3 and 5 the impact of k on flux values is modest and 
varies smoothly across a plausible range of k- values [0,2]. 
Nevertheless, its effect is measurable in light of the provided 
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formulas, contributing an extra tool to afterglow fitting and 
addressing the nature of GRB progenitors. 

Beyond the context of GRBs, the provided prescriptions 
are useful for modelling synchrotron emission from spherical 
adiabatic blast-waves of arbitrary velocity (with the limita- 
tions analysed in section 5.5) as long as they have swept 
up enough mass to be decelerating. Obvious applications 
include type Ibc supernovae (Soderberg et al. 2010), often 
associated with GRBs (Woosley & Bloom 2006) and mildly- 
or sub-relativistic spherical outflows from binary neutron 
star mergers. The latter are candidates for providing the 
EM counterpart (peaking at radio frequencies) to a possible 
signal of gravitational waves (Metzger & Berger 2012). By 
applying the ST scalings of the presented flux-prescriptions 
on mildly relativistic outflows we show that prospects of de- 
tecting such radio signals from within the horizon of gravita- 
tional wave detectors, LIGO (Abbott et al. 2009) and Virgo 
(Acernese et al. 2008), are realistic (Nakar & Piran 2011). 

It is interesting to note that apart from the dependence 
on p and k, we have also found that the sharpness of a spec- 
tral break can be influenced to some extent by the micro- 
physical parameters e c , es and £. This has been particularly 
seen in spectral breaks that involve absorption. The reason 
for this is the dependence of the absorption coefficient on the 
chosen microphysics through eq. (15). The microphysical pa- 
rameters in effect regulate the physical depth of the blast- 
wave corresponding to a given value of the optical depth. 
Therefore an increase/decrease of results in a less/more 
diverse sample of local electron distributions contributing to 
the flux across a spectral break and thus a sharper/smoother 
transition. We have chosen not to include the effect of the 
microphysics on the sharpness-formulas, as it typically influ- 
ences s by no more than 10 — 15% (the flux to a lesser extent) 
and it would greatly complicate the heuristic equations we 
present. 

Contrary to the approach on the evolution of a spectral 
break's sharpness, where the introduction of £nr is useful, 
when describing the temporal evolution of the spectrum's 
critical parameters we deliberately choose not to use such 
a timescale. The reason is that, as it turns out, there is no 
such thing as a single global timescale applicable to the be- 
haviour of all observable quantities. Instead, every critical 
parameter of the spectrum is chracterised by its own break 
time, the meeting point of the BM and ST asymptotes. One 
can verify that by computing to of eq. (27) for a few critical 
parameters of the same model, by equating the asymptotic 
expressions. They will be found to differ by factors up to 
a few. This happens because at any given observer time all 
these parameters are affected by contributions of radiation 
from various parts of the outflow, emitted within a range of 
lab-frame times. For each parameter the weight of these con- 
tributions will differ, leading to the inference of contrasting 
timescales by an observer. This stresses the need for models 
that can naturally account for this kind of features, by im- 
plementing accurate calculations of the blast- wave dynamics 
and the shape of the spectrum. 

By inspection of Fig. 4 one can realise that the bro- 
ken power-law approach is an approximation to the ac- 
tual behaviour of any critical spectrum parameter during 
the transrelativistic phase. The parameter that exhibits the 
strongest deviation from this description is v^- The reason 
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Figure 10. Simulation-based lightcurve taken at 10 20 Hz (no 
cooling taken into account). A bump at ~ 500 days is appar- 
ent. For comparison we have plotted the constructed lightcurve 
based on the presented flux prescriptions. Model parameters arc: 
E 52 = 1, n = 1, V = 2.3, k = 0, £ = 1, e e = 10- 4 , e B = 
10- 2 , d 2S = 1, z = 0.56. 

for this can be traced to the behaviour of the flux in the 
optically thin part of the spectrum. An example of this be- 
haviour is shown in Fig. 10, where a feature readily apparent 
is a smooth bump centered at ~ 500 days (this can also be 
seen in Fig. 10 of van Eerten et al. 2010a). This introduces 
a similar feature in the temporal profile of ^ a 2. As a result, 
the actual values of that frequency can deviate as much as 
~ 15% from the fitting function (smooth power-law break) 
at observer times relatively close to £nr- This can have an 
effect on constructed lightcurves if the self-absorbed part 
of spectrum 2 is used and only during the transrelativistic 
phase. The impact on flux levels is stronger than the devia- 
tions shown in Fig. 10 because the flux at the v 2 and u 5 ^ 2 
segments of spectrum 2 scales as i/~^ p+4 ^ 2 . We therefore 
recommend using eq. (33) in all cases as this method pro- 
vides a more accurate and consistent way of constructing 
spectra and lightcurves. 

By now there are a number of studies in the literature 
that present formulas calculating spectra from GRB after- 
glows. The importance of taking into account the blast-wave 
structure and the exact shape of the spectrum has been 
stressed by discrepancies in the derived values of physical 
parameters between simple (Wijers & Galama 1999) and 
more elaborate (Granot & Sari 2002; van Eerten & Wijers 
2009) models. Inclusion of details regarding the shock struc- 
ture can be done either analytically (in either of the two 
asymptotic regimes of the dynamics) or through the perfor- 
mance of simulations, as is done in the present work. The 
advantage of numerical simulations is that they cover the 
transrelativistic phase of the outflow and the level of detail 
they provide in all cases. The disadvantage is the price they 
come at, both in terms of time and resources. 

In this paper we provide an efficient way of utilizing 
the benefits of simulations, as those are reflected on the pre- 
sented analytic prescriptions. A similar approach has been 
taken by van Eerten & MacFadyen (2012), who are basing 
their fitting method on the scale-invariance of light-curves. 
This method naturally accounts for features like sideways 
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spreading of the jet and off-axis observation angles, features 
that only arise in simulations of at least two dimensions and 
cannot be captured in the context of this research. However, 
it requires the use of a large database of light curves which 
does not yet exist. So far no study of the afterglow radiation 
using 2D hydrodynamic simulations (Zhang & MacFadyen 
2009; Wygoda et al. 2011; De Colle et al. 2012) has resulted 
in the derivation of flux prescriptions. In fact this is the first 
time, even for the simple spherical case, that simulation- 
based flux prescriptions beyond the BM phase are presented. 
The box-fit method of van Eerten et al. (2012) does provide 
a fitting code but requires the use of a parallel computer 
network in order to fit data by iterating through a "box" of 
simulations. Therefore, analytic flux-prescriptions based on 
ID simulations, as the ones presented here, can be the base 
for comparisons with future work in that direction based 
on 2D simulations. Moreover, ID models are always rele- 
vant both at observer times before the jet break (when most 
parts of the outflow are causally disconnected) and at late 
times when the outflow is roughly spherical and allow for ac- 
curate calorimetry of jetted outflows after the jet-break but 
well before spherical symmetry has been reached, as long as 
the observer is not far off-axis (Wygoda et al. 2011). 

In all the simulations that we have performed to arrive 
at the presented flux prescriptions, the microphysical pa- 
rameters have been kept constant throughout the run and 
at every part of the outflow. This is by no means guaran- 
teed and therefore introduces an uncertainty in our results. 
An interesting topic for further study is the implementa- 
tion of evolving microphysical parameters and their effect 
on the flux prescriptions. Such an evolution is expected on 
theoretical grounds for some of those parameters (Granot 
et al. 2006). A qualitative study of the evolution of £ at 
the shock front and the evolution of ee downstream has 
already been presented in van Eerten et al. (2010a). Mean- 
while, there is also growing amount of observational evidence 
for this process taking place in GRB afterglows (Panaitescu 
et al. 2006; Kong et al. 2010; Filgas et al. 2011). Incorpo- 
rating effects like time-dependence of the microphysics into 
flux prescriptions can extend the predictions of the standard 
fireball model and thus broaden the theoretical framework 
within which observations are currently being interpreted. 



7 CONCLUSIONS 

We have used high-resolution ID hydrodynamic simulations 
to calibrate flux scalings of synchrotron, self-absorbed radi- 
ation for GRB afterglows in the relativistic and Newtonian 
dynamical phases (BM and ST, respectively). The transition 
from the former to the latter is well described by approxi- 
mating the evolution of spectral parameters (maximum flux 
and positions of critical frequencies) by power-law breaks 
connecting the two asymptotic behaviours. The properties 
of these breaks have been modelled in terms of the values 
of the physical parameters describing the blast-wave. This 
way we have managed to encapsulate the precision of the 
performed simulations into a set of analytic formulas that 
trace the full evolution of GRB afterglows, from the ultra- 
relativistic to the Newtonian phase. Due to the general na- 
ture of the prescriptions, they are applicable to any source 



characterised by emission of synchrotron radiation from an 
adiabatic blast-wave. 

A numerical code containing a practical implemen- 
tation of the results presented in this paper combined 
with a fitting code is freely available on request and on- 
line at http://www.astro.uva.nl/research/cosmics/gamma- 
ray-bursts /software / . 
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